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Abstract. Wc calculate a new gravitational wave 
background limit using timing residuals from PSRs 
J1713+0747, B1855+09, and B1937+21. The new limit is 
based on 17 years of continuous data pieced together from 
3 different observing projects: 2 at the Arecibo Observa- 
tory and 1 at the 140ft Green Bank Telescope. This project 
represents the earliest results from the 'Pulsar Timing Ar- 
ray' which will soon be able detect the stochastic back- 
ground from early massive black hole mergers. 



1. Introduction 

There are two categories of sources that could produce a 
detectable level of stochastic gravitational radiation from 
the distant Universe. One is of cosmological origin and 
the other involves galaxy evolution. First, low- frequency 
to relic gravitational radiation may be gener- 
ated during inflation, according to string theory models of 
the early Universe, via both early evolution of the extra 
dimensions and decay of cosmic strings. Cosmic strings, 
while now disfavored as the driver behind structure forma- 
tion in the Universe, still arise in many grand unified the- 
ories of particle physics after the epoch of inflation (Cald- 
well, Kamionkowski, & Wadley 1998; Hogan 2000; see also 
Maggiore 2000). Second, the coalescence of massive black 
holes (MBHs) during galaxy evolution could also produce 
a detectable level of Gravitation Waves (GWs). (See Lom- 
men & Backer 2001 for more on this.) Three cosmology 
experiments are being conducted or planned by a num- 
ber of investigators around the globe to either detect or 
place new limits on the stochastic background of gravi- 
tational radiation: polarization of the Cosmic Microwave 
Background Radiation (pCMBR), Pulsar Timing Array 
(PTA), and Laser Interferometer Space Array (LISA). 

The PTA consists of a regular observing schedule on a 
handful of the millisecond pulsars that are most accurate 
as clocks. So far, we have used the Arecibo Telescope in 
conjunction with the Green Bank 140ft telescope, but we 
plan on including the new Green Bank 100-m telescope 
that has just been commissioned and a collection of other 



European and Australian telescopes to maximize the ob- 
servation time and the range of observations available. 

Pulsar timing residuals are sensitive to gravitational 
radiation of wavelengths of about 1 yr, making the PTA 
complementary to LISA which will measure much shorter 
wavelengths, and the pCMBR which will measure much 
longer wavelengths. 

Using timing residuals from Arecibo Observatory ob- 
servations of PSRs B1937+21 and B1855+09 over 8 yr, 
Kaspi,Taylor, & Ryba (1994, hereafter KTR94) placed a 
limit of 



flgh'' < 6 X 10" 



(95%confidence) 



where fig is the fractional energy density in gravitational 
waves per logarithmic frequency interval and the Hubble 
constant Hq = lOO/i km s~^ Mpc~^. We have connected 
the KTR94 data to 10 years of Green Bank 140ft tele- 
scope data and 3 years of Arecibo Observatory data. The 
overlapping data sets span a total of nearly 17 years, dou- 
bling the baseline of Kaspi et al's results. The sensitivity 
of the PTA is proportional to 1 / {timespan)'^ , and there- 
fore we expect an increase in sensitivity of about 16 over 
the KTR94 results. 

First, in §||we describe the 3 different sets of observa- 
tions that have gone into calculating the limits presented 
in this article. In we connect the three data sets. In 
§1^ we show the calculation of the new gravitational wave 
limit. §^ compares the precision of pulsar clocks to those of 
atomic time standards. §0 is an update of the relationship 
between P and timing noise presented by Arzoumanian 
et al. (1994). In s|| wc discuss the possibility that a planet 
orbits pulsar PSR B1937+21. Finally in §|l^ we present 
our conclusions. 



2. Observations 

At the NRAOQ 140-foot (42.7 m) telescope we observed 
PSR J1713+0747, PSR B1855-H09, and PSR B1937-h21 in 



^ The National Radio Astronomy Observatory (NRAO) is 
operated by Associated Universities, Inc., under contract with 
the National Science Foundation. 
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addition to others 4-6 times per year at radio frequencies 
near 800 and 1400 MHz from 1989 October to 1999 July. 
For details on the analysis of these data see the report 
by Backer et al. (1993) on results from the first half of 
the data set. Observations of PSR J1713+0747 started 
after its discovery in 1994 (Camilo, Thorsett, & Kulkarni, 
1994). 

We have been conducting monthly observations at 0.43 
GHz, 1.4 GHz and 2.4 GHz of an array of MSPs using the 
NAIC0 Arecibo Observatory 300 m telescope since Decem- 
ber 1997. These data are used to make precision arrival 
time measurements for a variety of astrophysical goals. 
We used the Arecibo-Berkeley Pulsar Processor (ABPP), 
which is a multi-channel, coherent dispersion removal pro- 
cessor^ with 112-MHz total bandwidth capability (For a 
detailed technical description of the hardware see Backer 
et al.(1997.) For PSRs J1713-I-0747 and B1855+09 we use 
56-MHz bandwidth for observations at 1.4 GHz, and 112- 
MHz for 2.4 GHz, while for PSR B1937+21 we use 45- MHz 
bandwidth at 1.4 GHz and 56-MHz at 2.4 GHz. 

Calibrated total intensity profiles were formed from 
signals with orthogonal circular polarization. The profiles 
were then cross correlated with a template to measure 
times of arrival (TOAs) relative to the observatory atomic 
clock. Small errors in the observatory UTC clock, of order 
1 fis, were corrected based on comparison of local time 
to transmissions from the Global Positioning System of 
satellites (GPS). The templates used for cross correlations 
were constructed by fitting a set of Gaussian components 
to long-term averages of observations, using the software 
described in Kramer et al. (1994) and Kramer (1994). This 
model fitting scheme is described extensively in Lommen 
(2001). We use the model to generate noise-free templates 
with a specific common fiducial point at all frequencies. 
Lommen (2001) shows that the templates generated in 
this manner have at most 2 fis of error in the fiducial 
point between 1400 and 2380 MHz. 

In addition to the two data sets mentioned above, we 
also have incorporated the archival Arecibo data from 
KTR94 on PSR B1855+09 and B19374-21 into our anal- 
ysis. See the original paper for details on these data. 

3. Analysis 

As we just mentioned, TOAs are calculated calculated via 
cross-correlation with a template. Before these TOAs can 
be analyzed several corrections are required: 

TO A' = TO A + Aciock + Adm + A 

backend 

Aciock is comprised of a number of clock corrections which 
in our case transform a TOA which is referenced against 

^ The National Astronomy and Ionosphere Center Arecibo 
Observatory is operated by Cornell University under contract 
with the National Science Foundation. 

^ 'coherent' means that the dispersion is removed in the volt- 
age domain prior to power detection. 



an observatory maser, to one which is referenced against 
UTC. Adm causes TOA' to be the time that the pulse 
would have arrived had it been at infinite frequency (lower 
frequencies are delayed more by the ISM). Abackend cor- 
rects for any delays that are dependent on the observing 
backend used. In our case we have a digital latency which 
is dependent on channel bandwidth and a mid-scan cor- 
rection which is a result of our on-line folding. After all 
of these corrections, TOA' is independent of observatory 
clock, frequency, and observing backend. TOA' marks the 
space-time event of the arrival of a defined fiducial point 
of the pulse at the telescope on a defined time scale. 

4. Connection of Data Sets 

In the TEMPO analysis arrival time of pulses at the ob- 
servatory is referred to that at the barycenter of the solar 
system. The solar time corresponding to the rotation of 
the earth is known as UTl. UTl is tabulated by mea- 
suring the difference between UTl and UTC. UTC is a 
solar time scale that runs at the rate of TAI, interna- 
tional atomic time. BIPM tables of (UTl-UTC) are used 
to orient the observatory in inertial space and make the 
translation from observatory to earth-center. (UTl-UTC) 
is between 7 and 32 seconds over the course of our observa- 
tions, and is known to 0.1 ms at all times. This translates 
to an uncertainty in TOA of only '-^ 0.1 ns. To calculate the 
earth's position TEMPO uses standard JPL ephemerides 
DE405. 1950.2050. We know the geodetic position of each 
observatory to within meters. The delays that remain, that 
are most difficult to account for, arc the delays between 
the arrival of the pulse at the fiducial point of the antenna 
as defined by VLBI coordinates and the completion of its 
transmission through the 'back-end' of the receiving sys- 
tem. These are typically in the range of 100-3500 ns. At 
Arecibo, for example, we know that the travel time for a 
signal from the control room to the Gregorian receivers 
and back is approximately 7 /is Q largely resulting from 
the ^1200 ft of fiber-optic cabling through which the sig- 
nal must travel from the platform to the control room. 
This would give us a ~ 3.5 fis delay in the TOAs. Typ- 
ically, we ignore these delays because wc take data from 
a single observatory, with a single operating system, and 
any constant delay such as this is seemingly irrelevant. 

However, these delays become important for connect- 
ing multi-observatory data sets. If all such delays at all ob- 
servatories engaged in pulsar timing were known, it would 
be possible to have a "Universal Pulsar Timing Array 
(UPTA)" in which aU TOAs from all pulsar data from 
all over the world could be combined a priori in single 
meaningful data base. This would necessitate the need for 
uniform template fiducial point processing. For example, 
if the world's pulsar astronomers could piece together 17- 
year data sets on 10 different MSPs with comparable pre- 

* Mike Nolan, private communication, 2001 Aug 19 
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Table 1. Offsets, in /is, between data sets. 



Pulsar 


KTR94-GB=' 


GB-ABPP'' 


J1713+0747 




-3.2 


B1855+09 


-56.9 


1.4 


B1937+21 


-14.4 


-0.12 



GB=Green Bank data 
^ ABPP=ABPP data (post-upgrade Arecibo) 



cision, the sensitivity of the UPTA would go down by a 
factor of 1 /VTO- 

In our comparatively small experiment, we attempt 
this, but can only account for everything but the final 15- 
60 /xs, and so in the end we end up doing a bootstrap 
connection. We have 3 overlapping data sets. For each 
pulsar, the average separation between two overlapping 
sections is used as the offset between the two data sets. 
These offsets, which will be included in the ITOA format 
that we will make publicly available, are shown for each 
pulsar and for each adjacent data set in Table In the 
case of the Green Bank and ABPP data sets, the TOAs 
were generated using the same series of Gaussians with 
the same fiducial point, so all template discrepancies are 
removed to within 2 /xs (see Lommen 2001 for more on this 
subject). 2 ^s in fact, represents the maximum error in go- 
ing between 1400 and 2380 MHz. The data we have con- 
nected is centered from 1330 MHz to 1420 MHz, and the 
error roughly scales with the size of the bandwidth over 
which one must construct templates (see Lommen 2001). 
Therefore the maximum error in connecting our range of 
data, due to our template cross-correlation is roughly a 
tenth of 2 /j,s, or 0.2 ^s. In the case of the offset between 
the KTR94 data sets and the Green Bank data sets, we 
added, a priori, an amount of offset corresponding to the 
difference between our fiducial point and the fiducial point 
of KTR94, who used the peak of the template. These a pri- 
ori numbers were a phase of 0.443259 for B1855-I-09 and 
0.353460 for B1937+21. The numbers shown in the table, 
are the offsets in addition to these template offsets. The 
sign of the offsets are meaningful, e.g. the numbers under 
the column "GB-ABPP" are the quantities acquired from 
subtracting the average ABPP residual from the average 
Green Bank residual during the epoch while they overlap. 
The offsets between data sets that we processed ourselves 
(GB - ABPP) are all less than 8 ^s. 

In order to perform a weighted fit across multiple 
data sets, it is necessary to be very careful of the aver- 
age weighting of one data set relative to another. Said 
another way, with a single set of data, all taken with a 
single observing system and analyzed uniformly, one need 
not worry about, for example, overestimating all the er- 
ror bars by a factor of 3. However, when one is combining 
two such data sets, if one has overestimated error bars, 
and a second has underestimated error bars, the fit will 
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Fig. 1. (a) PSR J1713-h0747, (b) PSR B1855+09, (c) 
PSR B1937-I-21 and (d) PSR B1937+21 fit for P. The 
open circles arc KTR94 data, the filled circles are Green 
Bank data, and the open diamonds arc ABPP data. 



tremendously favor the latter. In order to account for this 
problem, we normalized the error bars in each data set by 
the total variance per day in those data, i.e. 



f^normalizcd — f^original ^ 



where L is the length of the data set, in days, and aoriginai 
are the error bars that were published by KTR or cal- 
culated by us, using the method created by Christophe 
Lange (private communication) which was adapted from 
Downs & Reichley (1983). 

We used the standard TEMPO package]^ to perform 
weighted fits to a, S, ^a, l-is, P, and P, in each of the 
3 pulsars. Additionally in J17134-0747 and B1855+09 we 
fit for the 5 Keplerian parameters and the Shapiro delay 
parameters, m2 and sini (See Stairs, 1998 for an excellent 
summary of the process of fitting for the Shapiro delay). 

We show the resulting residuals in Figure In (a) , (b) , 
and (c), wc show J1713+0747, B1855-h09, and B1937+21, 
respectively, fitted for the parameters not marked with 
footnote "c" in the table. Additionally, in panel (d), we 
show B1937-I-21 fitted for P{v). 



http://pulsar.princeton.edu/tempo 
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In the next three sections, we discuss various imphca- 
tions of the residuals resulting from these fits. First, we 
assume the residuals shown in Figure are the result of 
GWs and use them to place new limits on the GWB. Sec- 
ond, we assume the same residuals are the result of timing 
noise intrinsic to the pulsars, and determine whether these 
pulsars are more or less noisy than expected based on the 
characteristics of the population at large. Finally, we con- 
sider the possibility that the residuals in B1937-I-21 are 
the result of a planet orbiting that body. 

5. Limit on Gravitational Radiation 

These data, representing a continuous 17-y set, provide a 
new limit on the level of background gravitational radi- 
ation present throughout the galaxy, and ostensibly the 
universe. In short, we assume that all the observed fluc- 
tuations in the timing residuals are due to the stochas- 
tic background of gravitational radiation. This yields the 
maximum possible gravitational wave spectrum present. 

To determine the limit they place, we first need to 
measure the spectrum of fluctuations that we observe in 
the timing residuals. The process of computing an ob- 
served noise spectrum from irregularly sampled data has 
been given much consideration by previous authors (Groth 
1975; Deeter & Boynton 1982; Deeter 1984; Cordcs & 
Downs 1985; KTR94). These authors are particularly con- 
cerned with "red" spectrum where imperfect sampling 
leads to coupling of spectral estimates at low frequency. 
Deeter (1984) compares the various methods and deter- 
mines that the method of fitting the data to orthonormal 
polynomials produces the most meaningful results. We 
therefore use this method, as did both Stinebring et al. 
(1990) and KTR94 in previous generations of this experi- 
ment. 

Stinebring et al. (1990) and KTR94 start with a calcu- 
lation which shows that the energy density of background 
gravitational radiation per frequency octave would have a 
frequency dependency oc and they go about putting 
limits on that assumed spectrum. Wc do the same. How- 
ever, Rajagopal & Romani (1995), Phinncy (2001), and 
Jaffe & Backer (2002) recently show that the energy den- 
sity per frequency octave would be proportional to 
rather than . Either way, this is steep compared to the 
intrinsic spectrum observed in young pulsars, and would 
not change our results significantly. 

A series of polynomials, Pi{t) where i is of the set 
0, 1, 2, 3 of corresponding order, is generated orthogo- 
nally over the sampling of the timing residuals. Wc started 
with the following set: 

i 

Pj=0,l,2.3(0 = 51*' 
1=0 



is found which minimizes the quantity 

«=i y j=o J 

where r(tn) is the measured residual at time <:„. The first 
three Cj's are covariant with the fitting of the phase, pe- 
riod, and period derivative in the pulsar model. C3, how- 
ever, is a measure of the amplitude of the spectral den- 
sity of the variance in the timing residuals at a frequency 
corresponding to 1/T where T is the length of the data 
set. Sm =< C"! > is the corresponding estimate of the 
power (variance) contained near frequency 1/T. In order 
to sample the spectrum of Sm over multiple frequencies, 
we divide the data into smaller and smaller subsets in 
time, by factors of 2. We refer to the divisor as the fre- 
quency index. A frequency index of 2 implies the data set 
was divided in half, and that therefore we are measuring 
frequencies of 2/T, or periods of about 8.5 years. 

The results of the spectral measurements are shown in 
the solid line Figure | for PSRs J1713-f0747, B1855+09, 
and B1937+21, using the fits corresponding to parts a, b, 
and c of Figure |l|, i.e. fits through iy only. Though this 
method works well for irregularly sampled data, the value 
of the SmS arc still dependent on sampling. In order to 
gauge the extent of this dependence we have simulated 
data sets of various spectral indices, with the same sam- 
pling as the actual data. Randomly distributed Gaussian 
noise was transformed into the Fourier domain, multiplied 
by a function with spectral index 0,2,3, or 5, normalized 
to 1 ^s^ y at a frequency of 1 y~^, and transformed back 
to a time series. Once in the time domain, the data were 
sampled identically to the actual data. 10,000 such realiza- 
tions were created for each of the 4 spectral indices. The 
power spectra of each was measured using the method of 
orthonormal polynomials described above. The results of 
this analysis are shown in dotted lines in Figure ^ for each 
of the 4 spectral indices. 

The average value of the spectral estimator, Sm, with 
spectral index 5 represents the estimate of Sm in the pres- 
ence of a gravitational radiation, < Sm{^gh^) >, where 
rig is the energy density in GWs, p, expressed as a fraction 
of the closure density, pc of the universe: 

_ _p _ SttGp 
' Pc 3H^o 

where Ho is the Hubble constant. We generally express fig 
in terms of h, i.e. Qgh^ where Ho = 100 h km s~^ Mpc~^. 
We need to translate the normalizing amplitude, (1 ps^ y 
at a frequency of y~^) into a particular value of Q.gh?. 
We use the relationship derived from the definition of Q,g 
(from Stinebring et al. 1990) 



The process of orthogonalization is done via standard 
Gram Schmidt reduction. A linear combination of the p^s 
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Table 2. Observed and computed spectral densities. 



12 4 8 



1 2 4 a 16 1 2 4 8 16 



Frequency index, m 

Fig. 2. log 5m vs frequency index. Solid line shows mea- 
sured values for each pulsar. Dotted lines show simulated 
data for spectral index 0. Dashed lines show simulated 
data for spectral indices 2, 3, and 5. (a) PSR J1713-I-0747, 
(b) PSR B1855+09, and (c) PSR B1937+21. Frequencies 
corresponding to frequency index are m/L where L is the 
length of the data set: (a) 9.2 yr, (b) 15.6 yr, and (c) 16.8 
yr. 



where Pg{f) is spectral density, i.e. the same quantity we 
are attempting to measure with Sm- This means that val- 
ues of < Sm{flgh'^) > from data with 1 /is y at a frequency 
of need to be divided by 752 in order to represent the 
spectral density when ilgh^ = 10~^ which is the value we 
represent in Table |^. 

Figure ^ shows that the measured spectra of both 
J1713-h0747 and B1855+09 are both quite flat com- 
pared to the spectral index 5 model. We would expect 
J1713-I-0747 to be somewhat flatter than the others sim- 
ply by virtue of the plot representing a higher frequency 
region of the spectrum. PSR B1937-I-21 demonstrates a 
steep spectrum, comparable to the spectral index 5 model 
spectrum, but this may be due to other sources as we will 
discuss later (e.g. see 

In order to calculate the expected influence of gravita- 
tional waves on the spectral measurements, Sm, we needed 
to know the contribution from purely white noise. To this 
end, we simulated residuals that were only influenced by 
our known quantities of instrumental noise. We created 
white noise at the same sampling as that of the three 
pulsars we are considering. The white data from each ob- 
servatory were normally distributed random numbers nor- 
malized such that their RMS matched the median of in- 
strumental noise from that pulsar at that observatory. The 
value used for instrumental noise was the standard devi- 
ation of the mean of the timing residuals about the mean 
for a particular day, i.e. 



Pulsar 


m 








J1713+0747 


1 


348 


538 


3.21 




2 


6.20 


435 


0.108 




4 


22.4 


410 


3.52 X 10"^ 




8 


4.17 


420 


1.14 X 10"** 




16 


111 


278 


3.34 X 10"* 


B1855-f09 


1 


84.8 


137 


121 




2 


120 


182 


5.01 




4 


15.6 


210 


0.195 




8 


26.6 


204 


5.57 X 10"^ 




16 


54.7 


155 


2.51 X 10"" 


B1937+21 


1 


29379 


0.025 


164 




2 


203 


0.028 


9.87 




4 


11.8 


0.035 


0.277 




8 


7.93 


0.035 


8.39 X 10^^ 




16 


3.43 


0.034 


2.87 X 10"" 



stdev 



N 



GB=Green Bank data 

ABPP=ABPP data (post-upgrade Arecibo) 



where rt is a single residual of which there are N averaged 
to form a daily average and u is the mean of the residuals 
on that day. We created 10,000 random realizations for 
each of the 3 pulsars, and measured their spectrum using 
the same technique shown above. The average value of 
these data we call < Sm >w 

The values of Sm, < Sm >w, and < Sm >g for each 
value of the frequency index, m, are given in Table ^. 

In order to attempt to rigorously calculate an upper 
limit to background gravitational radiation present in our 
residuals, we duplicate the analysis of Thorsett & Dewey 
(1996, hereafter TD96). There is significant controversy 
around this method, which we address following the re- 
sults. TD96 create a statistic, Stati ~ ■mSm/[< Sm >w 
+ < Smi^gfi^) >g], which is normally distributed with m 
degrees of freedom. The statistic is essentially the mea- 
sured value of the spectral estimator divided by the pre- 
dicted value, based on Monte Carlo simulations. 

We use a Neyman-Pearson test to determine whether 
the distribution of Stati, for various values of flgh^ is 
significantly different than the corresponding statistic in 
the absence of gravitational radiation. The validity of the 
Neyman-Pearson test is complicated, but it has a simple 
construction. Consider a particular measurable, A, with 
probability distribution SI. Correspondingly, consider an- 
other measurable, B, with probability distribution S2. The 
Neyman-Pearson test essentially tests the extent to which 
the distributions overlap. If they have considerable over- 
lap, they are indistinguishable, if not, they aren't. The 
quantitative test of 'considerable' involves computing the 
line which delineates 5% of the area under the curve, repre- 
senting the least likely outcomes of measuring A. Consider 
Figure ^ which shows two probability distributions. The 
area to the left of the dotted line is 5% of the total area 
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8x10 



" 6x10 



4x10 



^ -7 
o 2x10 ' 




-7 



8x10 



« 6x10 



4x10 



^ —7 
o 2x10 ' 




Fig. 3. An example of two different probability density 
distributions. 5% of the area under the solid curve is to 
the left of the dotted line. 



of the solid curve. The mean value of the dotted curve, 
falls just to the left of the dotted line. Thus, we can say 
that these two distributions are different from each other 
at the 95% confidence level. 

One of the subtleties of the Ncyman-Pearson test is 
that one of the two hypothesis is considered the 'null' hy- 
pothesis, and this is the hypothesis one accepts unless the 
data compels one to choose otherwise. Thus, the Neyman- 
Pcarson test, by conjecture has a preferred answer, i.e. 
whatever you choose as the null hypothesis. Our null hy- 
pothesis will be that there is no gravitational radiation in 
the universe. This is of course incorrect. This is the first 
issue McHugh et al. (1996) have with this method. We are 
inclined to agree that this is not a good construction. 

TD96 used a likelihood ratio, which is an unnecessarily 
complicated way to construct the Neyman-Pearson test. 
We merely need to compare the two distributions that are 
shown in the numerator and the denominator of the ratio 
TD96 use. 



5*0 — Ilm=l,2,4,SXm 



so represents the null statistic, and is shown in the 
solid line of Figure ^ SI represents the statistic when a 
certain amount of gravitational radiation flgh^ is present. 
We have plotted SI in Figure ^ when its value is equal to 
the 5% limit that we calculate, namely, ilgh"^ < 2.8 x 10~^ 
at a 95% confidence level. 

TD96 calculated flg/i^ < 1.0 x 10"® at a 95% confi- 
dence level using only the KTR94 B 1855+09 data. Using 
the technique of TD96 our expanded set of B1855-I-09 data 
actually place a weaker limit than the KTR94 data alone, 
because of the large scatter in the Green Bank data points. 
The technique separates the data into smaller and smaller 



Fig. 4. The solid line shows the probability density of the 
null statistic, SO, while the dashed curve shows the prob- 
ability density of the SI statistic at the 95% confidence 
value. The dotted line is the same as in Figure || 

sections in order to measure shorter frequencies. The aver- 
age of these shorter frequency sections is therefore falsely 
enlarged in our calculations. Using only the m = 1 term of 
our data we have: flgh"^ < 9.8 x 10~^ at a 90% confidence 
level. 

McHugh et al. (1996) rejected the analysis of TD96 be- 
cause of their use of hypothesis testing rather than param- 
eter estimation. They completed a full Baycsian analysis 
of the same data which resulted in a much less restrictive 
limit, flgh"^ < 9.3 x 10~*, than the analysis of TD96 using 
the same data. 

These opaque calculations cover up the essential na- 
ture of the results. After KTR94 we can make a much 
more simple estimate of the limit these data place on the 
GWB. We use KTR94 equation 6 which relates the energy 
density, p, in gm cm~^ of a gravitational wave to its affect 
on TOA in /is. A, and a frequency, f, in Hz. 



P 



V416G J 



The largest amplitude sinusoid that one could conceivably 
fit to the B1855-t-09 data has A = 3 fxs and / = 1/17 yr = 
1.9 X 10"^ Hz. This gives p = 3.2 x 10"^* gm cm"^, or 



usmg 



P_ 

Pc 



SnGp 



= 2 X IQ^^h- 



which is more than an order of magnitude smaller 
than the limit found by KTR94. We assume Ho = 
100 h km s-i Mpc"^ 

As explained in Lommen & Backer (2001), a gravity 
wave propagating through the galaxy would have an effect 
both on the emitting site (the pulsar) and the receiving 
site (the Earth). Neglecting geometrical considerations for 
a moment, this could imply that the limit we have just de- 
rived will on average be a factor of 2 too high, i.e., that 
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the response that we see in the pulsar timing, since it re- 
sults from effects at both sites, is essentially doubled. The 
factor is not so simple, but depends on (a) the direction 
of propagation and polarization of the impinging gravita- 
tional wave, and (b) the geometrical relationship of the 
pulsars to the Earth. For a stochastic background of grav- 
itational waves, however, we can calculate the expected 
effect. If the GW at the emission and reception sites con- 
structively interfere then C3 will increase by a factor of 
2, and Sm will increase by a factor of 4. If emission and 
reception destructively interfere then Sm will be close to 
(it will certainly be non-zero due to irregular sampling) . 
The expectation value, i.e. the ensemble average over the 
47r directions of propagation, and over the wavelengths 
in question is a factor of 2, so the limit is actually half 
what we, or any of these other groups, measure. This of 
course leaves aside the issue of the variance in the GWB. 
The background value we measure using these data place 
a limit on the GWB in a specific place and epoch: in our 
galaxy and during the 2 decades in which this experiment 
took place. Jaffe & Backer (2002) arc working on a real- 
istic version of the GWB which takes these specifics into 
account. 

We have assumed up until now, that PSR B1855-f 09 
produces the best limit on the GWB since it demonstrates 
the most stable timing and a long baseline. However, what 
if there was a geometrical situation in which B1937-I-21 
was experiencing a large perturbation due to GWs while 
a null occurred at PSRs J1713+0747 and B1855+09? We 
first need to consider that the limited range of data on 
J1713-I-0747 was preventing our detection of any cubic 
therein. In Figure ^ we performed the same fit that we did 
for part (c) of Figure but using only the limited range of 
the J1713+0747 data. The obvious cubic is gone but there 
is a 'sawtooth' present at the level of 2 /xs that we do not 
see in the J1713-I-0747, so we conclude that we are seeing 
no comparable gravitational wave in J1713-I-0747. 

To attempt to find a geometrical situation which pro- 
duces an approximate null in J1713-I-0747 and B1855-f 09 
we inspected the range of possible geometrical multiplica- 
tive factors, (1 — 7)/2 over the 47r sphere of possible 
gravitational wave directions for each of the three pul- 
sars. See Lommen & Backer (2001) for an explanation 
of this factor. Figure ^ shows the factor for 2n values of 
RA while DEC=0° for each of the 3 pulsars. B1937H-21 
is shown by the solid line. J1713-I-0747 is the dotted line, 
and B1855-I-09 is the dashed line. The three pulsars are 
close enough together in the sky that (1 — 7)/2 is highly 
correlated over the sphere. We chose the DEC=0 line to 
show because it essentially produces the most favorable 
spot on the sphere for a B1937-I-21 gravity wave enhance- 
ment with a near-null in J1713-I-0747 and B1855-h09. At 
DEC=80° J1713-h0747 has (1 - 7)/2 = 0.005, B1855+09 
has (l-7)/2 = 0.05, while B1937-f 21 has (l-7)/2 = 0.12. 
That gives B1937-I-21 a 2.4x enhancement over B1855-t-09 
and a 24x enhancement over J1713-I-0747. The relative 
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Fig. 5. Residuals from B1937-I-21 with no v removed, 
fitted only to the range of the J1713-f 0747 data. 

strengths of J1713-f 0747 and B1855-t-09 can be exchanged 
somewhat as you can sec in Figure ^. We conclude that it 
would be very difficult to arrange for the 'deviant' residu- 
als we see in B1937-I-21 to be produced by gravitational ra- 
diation without producing a similar, but slightly smaller, 
- 10 effect, in either J1713+0747 or B1855-I-09. In 
addition, a gravitational wave disturbance, which, when 
reduced by a factor (1 — 7)/2 = 0.12 still yields a 30 /is 
residual deviation is unrealizable using any known grav- 
itational wave producer in the universe (see Lommen & 
Backer 2001). 

Finally, we consider whether emission and reception 
site disturbances could be destructively interfering in the 
case of B1855+09 and constructively interfering in the 
case of PSR B1937-F21. The stochastic GWB can be 
thought of as "crinkled" space-time. We approximate the 
crinkling as a sine- wave of equal amplitude but arbitrary 
phase near each of the 2 pulsars and also near the earth. 
We begin by computing the enhancement of the amplitude 
of a sine wave by adding two equal sine-waves together, 
with a variable displacement, (/>, between the two. This en- 
hancement we call E{4)). Thinking of B1937+21 and the 
earth as emitter and receiver, we may, for example, calcu- 
late the probability that £{(()) will be 5 or greater: 0.001. 
To obtain the probability that B1937-|-21's perturbation 
will be enhanced over Bl 855-1-09 's perturbation by a fac- 
tor of 25 (roughly the ratio of their peak deviations) or 
greater we compute the quotient: 



Q(A0) 



for A(/) = [0, 2tt]. The probability that Q(A(/)) > 25, 0.016, 
is the probability that the disturbance in B1937-I-21 is 25 
times that in B 18554-09. Thus, the residuals we observe in 
B1855-|~09 render it highly unlikely that the large ampli- 
tude quasi-sinusoid observed in B1937-f 21 is the result of 
a GW. Consideration of PSR J1713-f0747 would suppress 
the probability further but not by an equal factor. 
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Fig. 6. (1 - 7)/2 vs RA for DEC=0 for each of the 3 
pulsars. B1937+21 is shown by the sohd line. J1713+0747 
is the dotted line, and B1855+09 is the dashed line. 

6. Fractional Stability of Terrestrial Atomic Time 
Standards vs. Fractional Stability of Neutron 
StEir Rotations 

At the baselines considered in this paper, namely longer 
than 10 years, the fractional stability of neutron star ro- 
tations rivals that of terrestrial atomic time standards. To 
make the comparison quantitative we used the statistic 
(Tz proposed by Matsakis, Taylor, & Eubanks (1997) for 
describing pulsar and clock stabilities. The recipe for com- 
puting Oz given in Matsakis, Taylor, & Eubanks (1997) is 
very complete. We only give a bare outline of the compu- 
tation process here. We divide the timing residual or clock 
comparison, A"(i), into smaller and smaller subsections of 
time, starting with the full length, T, and going to T/2, 
T/4, T/8, etc. To each subset we fit the function 

X(t) = Co + Ci(t - to) + C2(< - to)' + C3(t - t^f. 

Oz is related to the average of C3 by the following formula 



< c? > 



1/2 



where the angular brackets denote the average, and r 
is the length of the subset of data. Figure |^ shows the 
fraction stabihty, vs. dataspan for [UTC-GPS], [TAI- 
PTB], [TAI-USNO], and [PSR-TAI] for PSRs B1855-h09, 
B1937-I-21. To compute the error we used the estimate 
provided by Matsakis, Taylor, & Eubanks (1997), which 
underestimates the error in the case of data with signifi- 
cant red noise, as is evident in the pulsar data. 

Figure shows that PSR B1855-f09 beats the [TAI- 
USNO] time scale at the longest time periods, and is es- 
sentially equally as accurate as [TAI-PTB] at timescales 
of 19 years. Since the fractional stability of B18554-09 im- 
proves steadily with increased dataspan, we expect that 
as additional data on B1855+09 is acquired, this pulsar 
will easily beat [TAI-PTB] at ^^25 year timescales unless 
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Fig. 7. Fractional Stability of 3 terrestrial timescales, and 
2 MSPs. 

TAI and PTB improve or the pulsar becomes less stable. 
Note that TAI is significantly better than USNO alone 
even though TAI is an aggregate that contains USNO, 
and similarly, for TT(BIPM) which contains TAI. It is 
also interesting to note that the PTA allows the possibil- 
ity of ignoring the earth clock altogether; one degree of 
freedom in the PTA can be used to solve for time, i.e., 
PSR vs PSR. 

7. Update on Timing Noise as a Function of 
Period Derivative 

Pulsar astronomers have long wondered whether the noise 
they see in timing residual plots, such as we show in Figure 
|l|is due to intrinsic instability in the rotation of the pulsar, 
or is an error in the model, such as a missing planet or un- 
modeled DM variations. Arzoumanian et al. (1994) pub- 
lished a compendium of measurements of "timing noise" 
in a number of slow and millisecond pulsars, and suggested 
that there is a relationship between timing noise and P, 
i.e., the larger the P the larger the timing noise. This 
implies first that the noise is intrinsic to the pulsar, and 
second that for each pulsar there is a minimum RMS in 
the residuals that no amount of accounting for systemat- 
ics will lessen. This has obvious implications for the PTA, 
and our effort to achieve sub-//s RMSs on a collection of 
MSPs. The data presented here along with recently pub- 
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Fig. 8. The timing noise parameter. Ag, as a function 
of log P . Values used are courtesy of Zaven Arzoumanian, 
except those shown in Table |[ and 1957+20 (Arzouma- 
nian, Fruchter, & Taylor, 1994), B1534+12 (Stairs et al, 
1998). 

lished results on MSPs allow us to update the picture of 
timing noise in MSPs as a function of period derivative. 
This was most recently done by Arzoumanian et al. (1994) 
who quantified the noisiness of a pulsar using the noise pa- 
rameter, Ag, defined by 

A8 = log(^i-p|i3^, 

where t = 10^ s. Arzoumanian et al. (1994) quote the best 
fit result to the scatter plot 

As = 6.6 + 0.6 log P 

which brought forth the rather depressing notion that wc 
would never achieve better than 1 p^s accuracy for the 
MSPs owing to the fact that P for the population is ~ 
10~^" or higher. 

We are happy to report that some of the MSPs used 
to create the low P end of the spectrum have been deter- 
mined to have lower Ag's than were presented in Arzou- 
manian et al. (1994). Wc present the new MSP data in 
Table |^ and the corresponding plot in Figure ^. 

The value used for PSR B1937+21 is the best fit to v 
letting all the parameters vary. Arzoumanian et al. (1994) 
were only able to place an upper limit on Ag for PSR 
B1855+09 and so are we, but ours is significantly lower, 
-6.9, down from -6. This comes from the best-fit value 
of V of —8 ± 4 X 10^^^, where the uncertainty quoted is 



that given by TEMPO. This is significantly smaller than 
the value previously obtained by KTR94, v = — 1.0 ± 
0.9 X 10^^^. The value for v that we use to calculate Ag 
is different from the value one obtains by only allowing v 
to vary. With only 1 degree of freedom, the value is not 
an appropriate upper limit. For J1713+0747 the fit refines 
nicely to give v = -2.6 ± 0.2 x 10"^^. 

PSR J0437-4715 was discovered by Johnston et al. 

(1993) . This 5.8-ms pulsar is very bright and has been 
shown to be very good for timing. The RMS of the resid- 
uals, folded and averaged at the binary orbital period, is 
only 35 ns. 

Since a third body was discovered around PSR 
B1257+12, a new upper limit on j> < -1.35 x 10"^^ has 
been published by Wolszczan et al. (2000). The value is 
an upper limit since Wolszczan suggests the distinct pos- 
sibility of a fourth body. While these new measurements 
do not change the value of Ag significantly from the value 
reported by Arzoumanian et al. (1994), it is now clear that 
the value represents an upper limit. 

PSR B1534+12 has an updated i> < 6 x lO"^^ (I.Stairs, 
private communication) which yields Ag < —5.4. 

PSR B1957+20 is possibly infiuenced by DM varia- 
tions, and also by an instable orbital system, as shown by 
Applegate & Shaham (1994) and Arzoumanian, Fruchter, 
& Taylor (1994). On the basis of the relationship between 
DM and RMS DM shown by Backer et al. (1993) we would 
expect PSR B1957+20, with a DM of 29, to have an RMS 
DM of about 0.0002 cm~'^pc. We can calculate the ex- 
pected influence of these DM variations on the estima- 
tion of V which was done on single-frequency data at 430 
MHz. A DM of 0.0002 cm^'^pc corresponds to a timing 
residual of 5 /iS. A 5 /iS trend over the course of the ob- 
servations performed by Arzoumanian, Fruchter, & Taylor 

(1994) for example, could produce a false i> of ~ 1 x 10^^^ 
given the 5-year observation length. This v corresponds 
Ag = —4.6, which is the value attributed to B1957+20 
by Arzoumanian et al. (1994). In other words, the timing 
noise measured by Arzoumanian et al. (1994) could be 
entirely due to unmodelcd DM variations. We therefore 
make this point an upper limit and emphasize that the 
true value is probably much smaller. 

It seems, then, that a significant fraction of MSPs fall 
below the line fitted by Arzoumanian et al. (1994). Either 
the relationship is not well described by a line, or the line 
has a steeper slope than was published. Either way, this 
bodes well for the PTA, which relics on submicrosccond 
accuracy in a handful of pulsars. 

Figure ^ is somewhat misleading in that it makes the v 
measured in B1937+21 look average. Many of the i>'s mea- 
sured are from single-frequency timing, so DM variations 
may be contributing to what we are calling 'timing noise.' 
There is, in fact, no MSP that displays such a convincing, 
significant v. We measure v in B1937+21 to 4 significant 
digits, where as most other MSP measurements are up- 
per limits or at best known to 20%. It is possible that 
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Table 3. Timing Noise Parameter vs v for the MSPs for which they are available. 



Pulsar 


P 


V (s-^) 


As 


Source 


J0437+4715 


5.72906(5) 


< 5 X 10"™ 


<-6.3 


vs'' 


J1012+5307 


1.7134(1) X 


-9.8(2.1) X 10"^'' 


-5.1 


la^ 


J1022+1001 


4.341(4) X 10"^" 


< 1 X 10"^' 


<-5.6 


kr" 


B1257+12 


11.4223(7) X 10"™ 


< -1.35 ±0.04 X 10~™ 


<-3.9 


wz^ 


J1713+0747 


8.54 X 10"''' 


< -2.6 ±0.1 X 10"^' " 


-5.6 


tw'^ 


B1855+09 


1.78 X 10~'** 


1 ± 6 X 10"™ 


< -6.9 


tw' 


B1937+21 


1.05 X lO"''-" 


1.515 ±0.001 X 10"™ 


-5.4 


tw' 


J2051-0827 


1.2737(5) X 10"™ 


< 2 X 10"™ 


<-3.5 


dos 



^ vs=van Straten et al. (2001) Although van Straten et al. (2001) did not calculate a limit i> we have estimated this quantity 
by taking the maximum deviation of their timing residuals (100 ns), converting this to a phase, and dividing by the cube of 
length of their data set, 3.4 y 

^ la^Lange et al. (2001) 

^ kr=Kramer et al. (1999) See note above''. We do the same for this pulsar using 50 ^ls and 1700 d. 
wz^Wolszczan et al. (2000) 

To calculate an upper limit we have allowed all the other parameters to vary along with i>. 
^ tw=this work 
^ do=Doroshenko et al. (2001) 



B1937±21 is displaying noise whose source is internal to 
the neutron star crust. Alternatively, perhaps B1937±21 
has a magnetic field structure that is evolving with a 20-y 
timescale, but we see such evolution at this magnitude in 
no other object. As an alternative we consider the possi- 
bility that B1937±21 is hosting a planet. 

8. Planet Around PSR B1937+21? 

The residuals shown in part (c) of Figure |l| are the result 
of a best-fit model a, 5, fia, iig, P, and P, and includes 
daily corrections for DM. One possible source of their cu- 
bic structure is the existence of a planet of binary period, 
Pb = 17.6 ± 0.2 y which yields an orbital separation of 8 
AU, using Kepler's Law with a central mass of 1.4 Mq. 
The residuals after including the planet in the model are 
shown in Figure ^. Asinj is 27.1 ± 0.7^s, yielding a com- 
panion mass of 0.08/ sini M0, where i is the inclination of 
the orbit to the line of sight. If real, this planet would be 
the smallest known extrasolar planet besides the lunar- 
mass planet around PSR B1257+12 (Wolszczan et al., 
2000). A planet of similar mass (0.15/ sini M^) has been 
found around PSR B1257+12 but is much closer to the 
pulsar (0.20 AU vs 8 AU) Wolszczan (1994). 

The planetary model adds 3 free parameters to the 
standard fit. To, Pb, and Asini. We compare the goodness 
of this model to adding the three parameters, i), v , and i' . 
The best-fit planetary model fits the data minutely better 
than the best i) through V fit (RMS of 1.14 ^s vs 1.18 
^s). The best fit parameters arc i> = 1.07(2) x 10~^^s~^, 
i>= 1.68(5) x 10-34s-4, and i> = -1.42(3) x lO^'^^s'^. We 
show in Figure |l^ the as a function of orbital period, 
and RMS as a function of orbital period, both of which 
show a minimum near 17 years, which is the length of our 
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Fig. 9. PSR B1937±21 residuals after fit for a planet. 

data set. The dotted line in both parts of the figure shows 
the value of the quantity when the fit to is turned on 
instead of the addition of the planet. Only i/, i>, Tq, and 
asini were allowed to vary in the fits. We did not allow 
a, 6, fia and /x/3 to vary for consistency with the 2-step 
fitting process described for B1937+21 in 

It is possible that the residuals you see in Figure |^ 
are the result of some intrinsic property of the pulsar (see 
^1^). For years this pulsar was thought to have significant 
"timing noise" causing the phase to wander as it does 
(KTR94). We were concerned that there may always be 
a well in the figure of vs Pb (e.g. Figure 10) where Pb 



is roughly equal to the length of the data set. We tested 
this possibility by doing the following. We simulated an 
extended data set, N years into the future, by assuming 
that i> is fixed. We then fit a planet to each data set and 
looked at the position of the well (in other words, the 
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Fig. 10. (a)^^ vs orbital period, aird (b) RMS vs orbital 
period for the model using a planet around B1937+21. 
The dotted line shows the value of the quantity when i) is 
used instead. 

best-fit Pb) vs N. What we found was that a 'planet' will 
produce a constant iy up through the year 2020, which is 
as far as we tested. There is a correlation between orbital 
period and the length of the data set, as shown in Figure 
[lT| . Some of the fits are not stable, such as the point at 
data span=18 y which yields a 3-y period. Notice that in 
general, the best-fit planetary orbit is much longer than 
the data set. This is easily understood if you consider that 
the planetary model needs to imitate the i) and can do 
so by extracting a fraction of an orbit from a putative 
planet. The non-continuous jumps to a higher Pf, shown 
in Figure |l^ are doublings of the orbital period. The 
and RMS for these simulated fits remains very similar to 
what we observe. We conclude that we could be fooled 
by a constant into supposing the presence of a planet. 
However, the previous section, §|^, shows that the large 
and highly significant i) in B1937+21 is strange relative to 
the population of MSPs. The only thing that will resolve 
this is additional ^ 10 years of data, assuming the orbit 
is about 20 years. 

PSR B1937-I-21 is regarded as having too large a P 
for being a recycled pulsar, i.e. its large P alone implies 
it is a 'young' pulsar although it is generally regarded 
as old Backer et al. (1993). We wondered if perhaps the 
unknown existence of an orbiting planet could have been 
skewing the measurement of P since the discovery of the 
pulsar. In fact, it could not. Our best fit P without the 
planet, and without assuming a i/ is slightly smaller, by 1 
part in ~ 10^ than our best fit P with the planet. 

9. Implications of Planet on Possible 

Evolutionary Scenarios of PSR B1937+21 

Millisecond pulsars are thought to be created in super- 
novae and spun up by accretion from a companion star 
during the red-giant phase of the companion. Forming 
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Fig. 11. Best- fit Planetary Orbital Period vs Length of 
Data Set for simulated extensions of PSR B1937-I-21 data, 
assuming a constant i). 

planets around such objects is complex and is the subject 
of much debate (Wolszczan 1998 and references therein). 
In fact, finding planets, such as this one, in similar mil- 
lisecond pulsars helps to solve the mystery in the creation 
of isolated millisecond pulsars. At the present, of the ~70 
known millisecond pulsars, 9 are solitary systems, meaning 
the presence of a companion has not been detected. Given 
that we suspect these objects have been spun-up by a com- 
panion, this presents an inconsistency. However, if planets 
are somehow the remains of an ablated or otherwise signif- 
icantly reduced companion, then we would expect to find 
similar planets around the other 8 known millisecond pul- 
sars. A planet with a 17-y orbit or longer, as we describe 
would not have been found in any other system, simply 
by lack of dataspan. 

10. Conclusion 

We have demonstrated that connection of data sets with 
sub-/is accuracy across multiple telescopes and many years 
is possible and yields high precision results for pulsar 
model parameters. The new limit that our data place on 
the energy density in background gravitational radiation, 
^ = 2 X 10~^/i~^, is below that yielded by the work of 
KTR94 by more than an order of magnitude. 

The data suggest the existence of a small (< 1 M0) 
planet around B1937+21. The best-fit orbital period is 
currently 17.6 y but that may change with the addition of 
more data. We have ruled out the possibility that the cubic 
present in B1937-|-21's residuals are due to a GW that has 
geometrically escaped detection in both J1713+0747 and 
B1855-I-09. We also show that if the cubic is due to a /> 
intrinsic to the pulsar, that B1937-I-21 is unique. We find 
this alternative therefore unlikely. 
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